%测试第一型曲线积分
clear; close all; clc;
f = @(x,y) x+y;
node = load("ex_node.mat");
ele = load("ex_ele.mat");
center = load("ex_center.mat");
node = node.ex_node;
ele = ele.ex_ele;
center = center.circucenter;
N = size(ele,1);
for i = 1 : N
    tri(i).nodes = [node(ele(i,1)+1,:);node(ele(i,2)+1,:);node(ele(i,3)+1,:)];
    tri(i).circucenter = center(i,:); 
    tri(i).edgelen = zeros(1,3);
    tri(i).edgelen(:) = compute_edge_length(tri(i).nodes);
    [tri(i).type, tri(i).which_angle] = judge_type(tri(i).edgelen);
    tri(i).normal = zeros(3,2);
    tri(i).normal(:) = get_normal(tri(i));
end
gauss_points = [0.3399810435848562648026658;0.8611363115940525752239465;
        -0.3399810435848562648026658;-0.8611363115940525752239465 ];
gauss_coef = [ 0.6521451548625461426269361;0.3478548451374538573730639;  
        0.6521451548625461426269361;0.3478548451374538573730639];
%node1 = [0,0];
%node2 = [2,2];
tmp_node = [node(ele(1,1)+1,:);node(ele(1,2)+1,:);node(ele(1,3)+1,:)];
node1 = tmp_node(1,:);
node2 = tmp_node(2,:);
vec = zeros(1,2);
vec(1) = node2(1) - node1(1);
vec(2) = node2(2) - node1(2);
judge = vec*tri(1).normal(1,:)';
k = (node2(2) - node1(2)) / (node2(1) - node1(1));
a = min(node1(1), node2(1));
b = max(node1(1), node2(1));
int1 = integral(@(x) (1 + k).*x - k*node1(1) + node1(2), a,b);
int1 = int1*sqrt(1 + k^2);
int = edge_int(f, node1, node2, gauss_points, gauss_coef);
x_point = [tmp_node(1,1),tmp_node(2,1),tmp_node(3,1)];
y_point = [tmp_node(1,2),tmp_node(2,2),tmp_node(3,2)];
x1_point = [center(1,1), center(1,1), center(1,1)];
y1_point = [center(1,2), center(1,2), center(1,2)];
plot(x_point, y_point,'-o','Color','b');
hold on;
quiver(x1_point',y1_point',tri(1).normal(:,1),tri(1).normal(:,2));
function int = edge_int(f, node1, node2, points, coef)
k = (node2(2) - node1(2)) / (node2(1) - node1(2));
h = abs(node1(1) - node2(1));
x_mid = 0.5*(node1(1) + node2(1));
y = @(x) k.*(x - node1(1)) + node1(2);
f_value = zeros(1,4);
for i = 1 : 4
    tmp_x = x_mid + 0.5*h*points(i);
    f_value(i) = f(tmp_x, y(tmp_x));
end
int = f_value*coef;
int = int*sqrt(1 + k^2)*0.5*h;
end

function [len12,len13,len23] = compute_edge_length(node)
len12 = compute_distance(node(1,:), node(2,:));
len13 = compute_distance(node(1,:), node(3,:));
len23 = compute_distance(node(2,:), node(3,:));
end

function len = compute_distance(node1, node2)
len_square = ((node1(1) - node2(1))^2) + ((node1(2) - node2(2))^2);
len = sqrt(len_square);
end

function [type, which_angle] = judge_type(len)
% 1 锐角， 2 直角， 3 钝角
% which_angle :1 12最长，2 13最长， 3 23最长
[max_len, which_angle] = max(len);
min_len = min(len);
len_sort = sort(len);
tmp_len = len_sort(2);
judge = min_len^2 + tmp_len^2 - max_len^2;
if judge > 0
    type = 1;
elseif abs(jugde) < 1e-13
    type = 2;
elseif judge < 0 
    type = 3;
end
end

function  normal = get_normal(tri)
normal = zeros(3,2);
index = 1;
for i = 1 : 2
    for j = 2 : 3
        if ( i ~= j)
            x_mid = 0.5*(tri.nodes(i,1) + tri.nodes(j,1));
            y_mid = 0.5*(tri.nodes(i,2) + tri.nodes(j,2));
            normal(index,1) = x_mid - tri.circucenter(1);
            normal(index,2) = y_mid - tri.circucenter(2);
            index = index + 1;
        end
    end
end
if tri.type == 1
    return;
elseif tri.type == 2
    max_edge_index = tri.which_angle;
    switch max_edge_index
        case 1
            node1 = tri.nodes(1,:);
            node2 = tri.nodes(2,:);
        case 2
            node1 = tri.nodes(1,:);
            node2 = tri.nodes(3,:);
        case 3
            node1 = tri.nodes(2,:);
            node2 = tri.nodes(3,:);
    end
    k = (node2(2) - node1(2)) / (node2(1) - node1(1));
    normal(max_edge_index,:) = [1,k];
    tmp = zeros(1,3);
    index = 1;
    for i = 1 : 2
        for j = 2 : 3
            if(i~=j)
                tmp(index) = normal(i,:)*normal(j,:)';
                index = index + 1;
            end
        end 
    end
    if min(tmp) > 0 
        normal(max_edge_index,:) = [-1,-k];
    end
elseif tri.type == 3
    max_edge_index = tri.which_angle;
    normal(max_edge_index,:) = -normal(max_edge_index,:);
end
end